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Abstract 

Employing extensive quantum-chemical calculations at the DFT/B3LYP and MP2 level, a quar- 
tic force field of isolated N-methylacetamide is constructed. Taking into account 24 vibrational 
degrees of freedom, the model is employed to perform numerically exact vibrational configura- 
tion interaction calculations of the vibrational energy relaxation of the amide I mode. It is found 
that the energy transfer pathways may sensitively depend on details of the theoretical description. 
Moreover, the exact reference calculations were used to study the applicability and accuracy of (i) 



O 

o 

Ph. 



> 

OO 

m 



^ \ the quasiclassical trajectory method, (ii) time-dependent second-order perturbation theory, and 
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(iii) the instantaneous normal mode description of frequency fluctuations. Based on the results, 
several strategies to describe vibrational energy relaxation in biomolecular systems are discussed. 



I. INTRODUCTION 

Amide modes in peptides and proteins are important probes to establish a structure- 
spectroscopy relation for biomolecules.-'^ The amide I modes, which are localized at the 
C=0 bond with frequencies of 1600—1700 cm~^, have attracted special attention because 
they are sensitive to the secondary structure of peptides.- This work is concerned with the 
vibrational energy relaxation of amide I modes. Recent experimental studies on model pep- 
tides as well as various small globular peptides have shown that the i/ = 1 ^ population 
decay time Ti of the amide I mode is about 1 ps for all systems considered.^A^ To explain the 
ultrafast energy and phase relaxation of amide I modes, a number of theoretical formulations 
have been giYQii,—^^'^^^^^^^^^^^^^'^^^^'^^'^^ many of which focus on N-methylacetamide 
(NMA, H3C-COND-CH3), a peptide-like small molecule containing only a single amide I 
mode. Employing standard biomolecular potential energy functions (such as CHARMM^i 
or GROMOS^^), quasiclassical trajectory simulationsi^ as well as time-dependent perturba- 
tion theory calculations^S, have been performed, qualitatively reproducing the subpicosecond 
relaxation time found for NMA in D20.-i^i^ 

Alternatively, various groups considered ab initio based potential energy functions to 
describe the dynamics of amide I vibrations.-'^'^i'i^iiiiii^'i^i^'i^ii^ Employing density functional 
theory (DFT) calculation at the B3LYP/6-31G+(d) level, we recently constructed a quartic 
force field of isolated NMA.— The dynamics of the amide I vibration was studied using 
the vibrational configuration interaction method,— '22i'2ii2^ including 24 vibrational degrees 
of freedom. The study indicated that the energy transfer pathways may sensitively depend 
on details of the theoretical description. Hence, one goal of the present paper is to examine 
the sensitivity of amide I relaxation with respect to the ab initio level (DFT/B3LYP and 
MP2), basis set (6-31G+(d) and 6-31G++(d,p)), and description of the solvent (none or 
continuum model). 

As the vibrational configuration interaction method gives a numerically exact quantum- 
mechanical description of the vibrational dynamics of the model system, it provides an ideal 
means to test several approximations usually employed in the description of vibrational en- 
ergy relaxation.— The focus of this work is therefore to study the applicability and accuracy 
of (i) the quasiclassical trajectory method,— (ii) (time-dependent) second-order perturbation 
theory,— i^I and (iii) the instantaneous normal mode description^^ of frequency fluctuations. 
Based on the results, we discuss several strategies to describe vibrational energy relaxation 



in biomolecular systems. 



II. THEORY AND METHODS 



A. Quartic force field 



While there are various ways to construct ab initio potential energy surfaces (PES),— ^^S 
we employ here a simple but accurate method called quartic force field, which is useful 
if we consider vibrational dynamics around a single equilibrium state of a molecule.— i^i^ 
In this approach, the PES is approximated by a fourth-order Taylor expansion around the 
equilibrium geometry 

V{{Qi}) = 7: /_^^kQk + TTy /_^ tklmQkQlQm + Tj ^_^ UklmnQkQlQmQn, (1) 

k k,l,m k,l,m,n 

where Qk and ujk denote the fcth normal coordinate and the associated harmonic frequency, 
respectively, and the coefficients t^im and Ukimn represent the third- and fourth-order deriva- 
tives of the PES. The above expression can be recast in the form of the n-mode coupling 
representation (nMR)^'2ii2^ 

V{{Q^}) = V'''^{{Q,}) + V'''^{{Q,}) + V'^'^iiQ,}), (2) 



V'^'^iiQ.}) = >J i-u^Q', + -t,,,Qf + ^uuQt) , (3) 

v^2Afi?,(|g^l) = Y. {It^nQ^Q' + l^^jnQ^Q') + E l^^v^Q'Ql (4) 

V'^'''{{Q^}) = Yl U,kQ^Q,Qk + Y. l^mQlQjQ'^- (5) 

i<j<k ij^j<k 

Note that we have omitted the 4MR terms V*^'^^{{Qi}) = YliiKjKkKi it'^ijkiQiQjQkQi- This 
strategy was successfully applied to several molecules.— i^i^i^ By adding the normal mode 
kinetic energy K = ^^ -P^ /2, we obtain the full approximate vibrational Hamiltonian H = 
K + V{{Qi}) for the system. 

In this study, the quartic force field was derived using the B3LYP or MP2 methods 
with the 6-31G+(d) or 6-31G++(d,p) basis set. First, the equilibrium structure and the 
harmonic frequencies were obtained for trans-NMA using the GaussianOS program package.— 
Then, the third and fourth-order derivatives were calculated by numerical differentiation of 
the analytic HessiaiiS^ using the SiNDO code developed by Yagi^ or by the "anharmonic" 
option of GaussianOS. 



B. Vibrational configuration interaction method 

In a first step, we calculated the vibrational energies and eigenstates of the above de- 
scribed model of NMA using the vibrational configuration interaction (VCl) method de- 
veloped by Bowman and coworkers.— To generate basis functions for the VCI calculations, 
we first performed a vibrational self-consistent field (VSCF) calculation of the vibrational 
ground state. The VSCF wave function is expressed as a direct product of one-mode func- 
tions or modals^ as 

*r"^=n€wo, (6) 

where n and / denote the vibrational quantum numbers and the number of degrees of 
freedom, respectively. The modals are determined by the VSCF equation 
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which yields the vibrational ground state (n=0) as well as the virtual modals. The VCI 
wave function is expressed as a linear combination of VSCF configurations 

m 

The VCI wave functions and energy levels are obtained by diagonalization of the VCI matrix 

All VSCF/VCI calculations were carried out using the SiNDO code^ for non-rotating 
molecules. The modals were expanded in terms of harmonic oscillator wave functions. The 
number of wave functions employed were 11, 9, 7, and 5 for {(j)^'^\ 0^^''}, {0^^-'- 0'-^^''}, {0^^^-*- 
0*^^^^}, and {0(^^)-0(^'')}, respectively. The mode index was labeled in increasing order of the 
frequency and the six lowest-lying modes were kept frozen. The VSCF configurations were 
selected by increasing energy until a cut-off energy of ~ 5000 cm~^, resulting in A'^ci — 6000. 
Once all eigenvalues {-En} and eigenf unctions {^&n^^} are obtained, it is straightforward 
to calculate the time-dependent wave function 

\m{t)) = ^(^;(^^i^(o))e-'^"*/'*i^;(c^), (10) 



which represents a numerically exact description of the above described quartic force field 
model. 

In all calculations discussed below we assume that the initial state is given by a single 
VSCF configuration \E'(0) = $y^^^, representing the n = 1 state of the amide I vibration. 
We thus obtain for the wave function 

with 

Oj(t) = $^Cj„aie-^"*/^ (12) 

n 

and -Pj(t) = |Oj(t)p for the time-dependent probability of state ^Y^^^. This allows us to 
define the harmonic energy of vibrational mode i as 

E,(i) = ^^.$^^jW^f, (13) 

j 

where the zero point energy of the mode was disregarded. Equation (fT3|) will be used below 

to discuss the vibrational energy relaxation in NMA. 

C. Time-dependent perturbation theory 

Employing time-dependent perturbation theory, recently Fujisaki et al. have extended 
Fermi's Golden Rule to the non-Markovian regime.—'^ As a stringent test of this formu- 
lation, we wish to compare the perturbative results to the exact VCI results obtained for 
NMA. Assuming zero temperature, the perturbative expression for the amide 1 ground state 
population is given by^ 

Po{t) = poo{t) ^ ^ > 7Z r^ 

- Y^ ^L/3 [1 - C0S{US -UJo.-UJ(3)t], (14) 

where ujs is the is the anharmonicity-corrected system frequency, Ua denote the (harmonic) 
frequency of the bath modes, and tsa/s represents the third-order coupling elements.—. More- 
over, we have introduced a Fermi resonance parameter— 
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which is proportional to 

V 



AE ■ '■''^ 

where \i) and |/) are the initial and final harmonic states, AV = V — Y2i^ uj1Q\/2 denotes the 
anharmonic vibrational coupling, and AE is the energy difference between \i) and |/). The 
Fermi resonance parameter is a key ingredient in the interpretation of vibrational energy 
transfer processes.— This concept has been useful for the classical dynamics characterization 
of the vibrational relaxation pathways in water— and in myoglobin.— 

D. Classical description 

Quasiclassical trajectory calculations are a well-established approach to approximately 
calculate the energy redistribution in molecular systems.— In this method, the initial state 
of the quantum system (e.g., the n = 1 state prepared by an infrared laser pulse) is rep- 
resented by a phase-space distribution (e.g., the Wigner distribution), which is sampled by 
an ensemble of classical trajectories. Starting by construction with the correct initial state 
and recalling that classical mechanics can be considered as a short-time approximation to 
quantum mechanics, a classical simulation should be a good approximation to quantum- 
mechanics in the case of ultrafast relaxation dynamics. ^^ 

In this work, we solved Newton's equation on the PES defined in Eq. ([T]), using Yoshida's 
sixth-order symplectic integrator with a 0.5 fs time step.— To generate classical initial con- 
ditions for the positions and momenta, we represent the normal modes {Qk-, Pk} in terms of 
classical action-angle variables {ri/j, (pk]^ 

Qk = \/{2nk + 'j)h/ujksm(f)k, (17) 



Pk = >y{2nk + j)hujk cos 4>k, (18) 

where the factor 7 = 1 accounts for the zero-point energy of the mode. To obtain the 
initial positions and momenta of the initially excited amide I mode, we associate the action 
Uk with the initial quantum state of the amide I mode, i.e., n^ = 1 for the first excited 
state. The remaining modes are for T = K initially in the ground state n^ = and 
only vibrate according to their zero-point energy. In all cases, the vibrational phases (pk are 
picked randomly from the interval [0,27r]. Using this initial conditions, the classical energy 



content of vibrational mode i is defined as 

E.{t) = ^^+'^{Qm'Er^ (19) 

where (...) represents an ensemble average over 100 nonequilibrium trajectories and Ef^^ 
denotes the zero-point energy of mode i. 

Although Eqs. (TT7|) and (fTSll represents a correct quasiclassical representation of the 
initial state of the system (in the sense that the vibrational energy distribution is correct), 
it may give rise to unphysical behavior due to the well-known zero-point energy problem of 
classical mechanics.— In quantum mechanics, each oscillator mode must hold an amount of 
energy that is larger or equal to the zero-point energy of this mode. In a classical trajectory 
calculation, on the other hand, energy can flow among the modes without this restriction. In 
the case of our ab initio model of NMA, for example, we have found that the high zero-point 
energy contained in the C-H stretch modes (~ 1600 cm~^) may be transferred to a large 
amount to low-frequency modes of NMA, thus yielding unphysical and unstable trajectories. 

Various approaches have been proposed to fix the zero-point energy problem,— however, 
most of these techniques share the problem that they manipulate individual trajectories, 
whereas the conservation of zero-point energy should correspond to a virtue^of the ensemble 
average of trajectories. Here, we adopt the method introduced in Ref. 
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which invokes 



quantum corrections to the classical calculation in order to restrict the classically accessible 
phase space according to the rules of quantum mechanics. At the simplest level of the 
theory, these corrections have been shown to correspond to including only a fraction 7 (0 < 
7 < 1) of the full zero-point energy into the classical calculation. The quantum correction 
7 can be determined by requiring that some observables of the system are reproduced,— 
e.g., that the amide I energy remains larger than the zero-point energy for all times under 
consideration. In the present case, we assigned full zero-point energy to the initially excited 
amide I mode, while for the remaining modes we choose 7 according to the classical thermal 
energy at T = 300 K. This protocol is simple to implement in standard MD codes and was 
shown to give good results in previous quasiclassical trajectory studies of vibrational energy 
redist ribution .— 



III. RESULTS AND DISCUSSION 

A. Quantum dynamics at various levels of ab initio theory 

The theoretical description of the amide I vibration of NMA first of all depends on the 
quantum-chemical parameterization of the quartic force field. Here, the following four levels 
of theory were employed: DFT/B3LYP with 6-31+G(d) or 6-31++G(d,p) basis set and MP2 
with 6-31+G(d) or 6-31++G(d,p) basis set. To facilitate the comparison to experiment, we 
first consider the vibrational frequencies of NMA, as obtained from the diagonalization of the 
VCI matrix i^. Figure [T] compares these frequencies to experimental values obtained under 
argon matrix conditions, ^^ focusing on vibrational modes between 1000 and 2000 cm~^. 
MP2 somewhat overestimates the anharmonic frequencies in comparison to DFT/B3LYP. 
Applying the 6-31+G(d) basis set, we obtain an overall deviation (mean absolute deviation) 
of 52 and 26 cm^^ at the MP2 and DFT/B3LYP level respectively. As the improvement 
due to the larger 6-31++G(d,p) basis is relatively small (compared to the increased effort), 
in the following we only show results using the 6-3H-G(d) basis. We confirmed that the 
energy transfer pathways are nearly the same in both cases. 
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FIG. 1: Vibrational frequencies of isolated N-methylacetamide (NMA), calculated from the quartic 
force field using the VSCF/VCI method. Compared are experimental^ and calculated frequencies 
at various levels of theory. The amide I vibration corresponds to mode ^ 23 and amide II to # 22. 



While all levels of theory at least qualitatively reproduce the vibrational frequencies of 
NMA, the deviations may be more significant when the vibrational energy relaxation of the 
amide I vibration is considered. Figure [2] shows the time evolution of the energy content [Eq. 
flT3|) ] of the initially excited amide I mode (# 23) as well as of the remaining modes of the 
molecule. As discussed previously,— in the B3LYP/6-31G+(d) calculation most of the amide 
energy goes to the 9th vibrational mode at 869 cm^^ via the pathway |23i) — » |92) and to the 
7th and 12th modes via the pathway |23i) -^ |7il2i). As discussed in Ref. Il6|, these normal 



modes correspond to motion of the C-terminal methyl-group atoms. The pathways can be 
explained by considering the Fermi resonance parameter Fsap defined in Eq. (TT5|) , which de- 
scribes within time-dependent perturbation theory the intensity of transition l^i) —>■ |ai/5i). 
Figure [3] shows the Fermi resonance parameters for all mode combinations (a,/5) together 
with its resonance factor |l/(ci;5 — uja~ uJjs)] and the anharmonic coefficient Its^/?!- Interest- 
ingly, it is the resonance factor rather than the anharmonic coefficient that determines Fsai3 
and thus the relaxation pathway. This underlines the importance of an accurate calculation 
of vibrational frequencies. 
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FIG. 2: Time evolution of the energy content of the initially excited amide I mode as well as 
the remaining modes of the NMA. Compared are quantum (red lines) and classical (green lines) 
calculations, obtained at the DFT/B3LYP (left) and MP2 (right) level of theory. 

Employing the force field obtained from the MP2/6-31G+(d) calculations. Fig. [2] shows 
a qualitative change of the amide I relaxation dynamics as compared to the results at the 
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FIG. 3: Fermi resonance parameter FsajS [Eq. ()15|) ] together with its resonance factor \\/{LOs~'-^a — 
u}pi)\ (cm) and the anharmonic coefficient |t5o/3| (kcal/mol/ A) . Compared are quantum calculations, 
obtained at the DFT/B3LYP (left) and MP2 (right) level of theory. 
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DFT/B3LYP level. While the same vibrations are involved in the process (i.e., mainly modes 
# 9, and # 7 and # 12), the amide I energy decreases only by ~ 20 % instead of ^ 80 %. 
Figure [3] reveals that this is caused by the fact that the value of the corresponding Fermi 
resonance parameters for MP2 are only half of the value of the DFT parameters. (Note that 
FsajS enters quadratically in the relaxation rate, which explains the overall factor of four.) 

As an approximative way to include solvent effects into the parameterization of the quartic 
force field, we have employed the implicit solvent model as implemented in GaussianOS— and 
repeated the above calculations. As an example. Fig. H] shows results obtained for implicit 
water and implicit methanol at the DFT/B3LYP level of theory. In both cases, one obtains 
a relatively weak amide I relaxation, which mainly involves the modes #7, #8, #11, and 
7^12. However, it should be stressed that in experimental^ as well as in simulations using 
explicit water solvenir^i^ there is a strong flow of vibrational energy into the solvent degrees 
of freedom which is, of course, not captured by an implicit model. Hence one should regard 
these calculations as another example of energy transfer dynamics on different PESs. 
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FIG. 4: Vibrational energy relaxation of NMA, obtained from DFT/B3LYP calculations using an 
implicit solvent model. Compared are quantum (red lines) and classical (green lines) calculations, 
obtained for implicit water (left) and implicit methanol (right). 
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B. Classical relaxation dynamics 

The numerically exact quantum calculations shown in Figs. [2] and H] provide a stringent 
test for quasiclassical trajectory simulations of vibrational energy transfer. The comparison 
of both methods in these figures reveals that for all cases considered the quasiclassical 
approximation yields at least qualitative agreement with the reference calculation. It is 
interesting to note, though, that the classical results deviate already at short times. This 
is a consequence of the fact that we haven't employed the full zero point energy in initial 
conditions fITTj) and (fTSi) . Doing so, we indeed obtain the correct dynamics for the first few 
hundreds of femtoseconds, but also observe unphysical flow of zero point energy at longer 
times (data not shown). ^- By employing full zero point energy in the initially excited amide 
I mode and thermal energy in all remaining vibrational modes, on the other hand, we obtain 
a reasonable classical approximation. 

Nonetheless, certain care should be taken when a quasiclassical approximation to quan- 
tum dynamics is employed. In classical mechanics, energy can flow among the vibrational 
modes without the restrictions of quantum mechanics. For example, one may observe classi- 
cal energy transfer, if the initially excited vibrational mode ujs is about half of the frequency 
of a high-frequency bath mode Ua (i.e., 2ujs ~ cu^), even though the initially excited mode 
contains only the energy of a single quantum. For the present force field, this behavior 
occurred when the amide II mode was excited and energy transfer to the C-H stretch modes 
was observed (data not shown). As in the case of the zero point energy problem, it is there- 
fore always advisable to check the main energy pathways of a classical calculations, in order 
to make sure that the classical energy flow could also happen in quantum mechanics. 

C. Perturbative description of relaxation dynamics 

Apart from classical approximations. Golden Rule-type methods as described in Sec. Ill CI 
have be proven useful to describe vibrational relaxation djTiamics.— i^ The comparison of 
this method to the reference calculations in Fig. [5] reveals that the perturbative expression 
gives quantitative results for the first few hundreds of femtoseconds. This may be expected, 
because the initial process is given by the energy transfer between the amide I mode and the 
doorway states (or first tiers) ,-^ which is correctly taken into account by the perturbative 
formula. At longer times, higher-order interactions and the subsequent relaxation of the 
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doorway states into other degrees of freedom may become important, which is not included in 
the second-order perturbative description. Nevertheless, in the present case the perturbative 
calculation yields qualitatively correct results for the first few picoseconds. 



T 1 1 r 

Full Quantum — 
Perturbation — 





FIG. 5: Comparison of exact and perturbative calculation of vibrational energy relaxation of the 
amide I mode of NMA, employing (left) B3LYP and (right) MP2 level of theory. 

A possible approach to extend the application of Golden Rule-type methods to the con- 
densed phase is to take into account the fluctuations of the solvent by calculating the in- 
stantaneous normal mode frequencies of the solute;^ To study the validity of this approach, 
we have calculated instantaneous normal mode frequencies for the present model, assuming 
both equilibrium and nonequilibrium initial conditions. As shown in Fig. [6], the instanta- 
neous normal mode frequencies of the amide I mode and other high-frequency modes of 
NMA may undergo fluctuations of several hundred wavenumbers. Even under equilibrium 
conditions, the amide I mode may vary by more than two hundred wavenumbers. This 
behavior is quite unrealistic, when we recall that, in experiment, typically fluctuations of 
tens of wavenumbers are observed. ^'^'^ We conclude that instantaneous normal mode fre- 
quencies cannot be regarded as "true" vibrational frequencies, when strongly anharmonic 
ab initio PES are employed. A more realistic description of the time-dependent frequen- 
cies is obtained, if we perform a geometry optimization of the solute molecule, leading to 
"quenched normal modes".— i^ Due to the absence of a fluctuating solvent, however, in the 
present case this procedure only leads to the global minimum of the molecule and therefore 
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to time-independent (standard) normal mode frequencies. 
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FIG. 6: Instantaneous normal mode frequencies of selected high-frequency modes of NMA, assum- 
ing (left) nonequilibrium and (right) equilibrium initial conditions, respectively. 

IV. CONCLUSIONS 

Following the previous work,— we have constructed various quartic force fields of iso- 
lated N-methylacetamide (NMA) and performed numerically exact VSCF/VCI calculations 
of the vibrational energy relaxation of the amide I mode of NMA. All levels of ab initio 
theory considered (DFT/B3LYP and MP2 with 6-31+G(d) or 6-31++G(d,p) basis set) re- 
produced the vibrational frequencies of NMA at least qualitatively. However, the various 
theoretical levels were found to lead to significant deviations of the vibrational energy relax- 
ation. This behavior could be explained through the Fermi resonance parameter [Eq. OlSp ], 
which includes the anharmonic coupling as well as the resonance condition of the energy 
transfer. The latter was found to predominantly determine the efficiency of the process. It 
should be noted that this sensitivity to the accurate description of vibrational frequency 
is particularly high in the case of isolated NMA. Larger peptides and proteins in aqueous 
solution provide a much higher vibrational level density, which facilitates efficient relaxation, 
since there are always some approximately resonant modes available.— Combined with the 
overall averaging effect of solvent-induced frequency fluctuations, this may render a more 
qualitative quantum-chemical description sufficient. 
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Taking the VSCF/VCI results as reference calculations, we have studied the validity and 
accuracy of various approximations usually employed. It was found that second-oder pertur- 
bation theory provides a faithful description of the first few hundreds of femtoseconds of the 
relaxation dynamics. This is because the initial energy transfer between the amide I mode 
and the doorway states is correctly described by the perturbative formula. At longer times, 
this approximation yields only qualitatively correct results, since higher-order interactions 
and the subsequent relaxation of the doorway states into other degrees of freedom may 
become important. To extend this method to the condensed phase, instantaneous normal 
modes are often invoked in order to account for the solvent-induced frequency fluctuations of 
the system. In the present case of strongly anharmonic ab initio PES, however, the resulting 
instantaneous normal mode frequencies were found to undergo unrealistically high fluctua- 
tions of several hundred wavenumbers. A more realistic description of the time-dependent 
frequencies may be obtained, if we perform a geometry optimization of the solute molecule, 
leading to "quenched normal modes" .—^ 

We have also performed quasiclassical trajectory simulations of vibrational energy trans- 
fer, which were found to be in at least qualitative agreement with the reference calculation. 
However, in classical mechanics, energy can flow among the vibrational modes without the 
restrictions of quantum mechanics. A well-known example is the zero-point energy problem, 
which in the present case was circumvented by employing full zero point energy only in the 
initially excited amide I mode and thermal energy in all remaining vibrational modes. More- 
over, in classical mechanics energy transfer may also occur to high-frequency bath modes. 
To assure that the classical energy flow also would happen in quantum mechanics, it is 
therefore always advisable to check the main energy pathways of a classical calculation. 
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